/****************************************************************************/
/* This file is part of FreeFem++.                                          */
/*                                                                          */
/* FreeFem++ is free software: you can redistribute it and/or modify        */
/* it under the terms of the GNU Lesser General Public License as           */
/* published by the Free Software Foundation, either version 3 of           */
/* the License, or (at your option) any later version.                      */
/*                                                                          */
/* FreeFem++ is distributed in the hope that it will be useful,             */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
/* GNU Lesser General Public License for more details.                      */
/*                                                                          */
/* You should have received a copy of the GNU Lesser General Public License */
/* along with FreeFem++. If not, see <http://www.gnu.org/licenses/>.        */
/****************************************************************************/
// SUMMARY : ...
// LICENSE : LGPLv3
// ORG     : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
// AUTHORS : Sala Lorenzo
// E-MAIL  : salallo80@gmail.com

// *INDENT-OFF* //
//ff-c++-LIBRARY-dep:
//ff-c++-cpp-dep:
// *INDENT-ON* //

#include "mode_open.hpp"
#include <iostream>
#include <cfloat>
#include <cmath>
#include <iterator>
using namespace std;
#include "ff++.hpp"
using namespace Fem2D;

/*!The class DxWriter permits to save in opendx format a "field"
 * (in the dx-language a "field" means the values of a function f(x,y,z) on a grid),
 * a time series (an ordered collection of "fields", so we have field0=f(x,y,z,t0), field1=f(x,y,z,t1)
 * and so on).
 * DxWriter creates two files: one with extension .data where it puts the position of the grid, the connessions,
 * the values; and one with extension.dx where it puts the time series.
 * Now you can save only scalar fields.
 * An example <code>
 * load "DxWriter"
 * mesh Th=square(5,5);
 * DxWriter ff("pippo");
 * Dxaddmesh(ff, Th);
 * Dxaddtimeseries(ff, "Vx",Th);
 * fespace Vh(Th, P1);
 * Vh vx=x*y;
 * Dxaddsol2ts(ff,"Vx",1.0, vx);
 * vx=0.5*x*y^2+0.2;
 * Dxaddsol2ts(ff,"Vx",2.0, vx);
 * cout<<"Ciao";
 * </code>
 */

class DxWriter {
	struct tsinfo {
		int imesh;	//!< index of the mesh
		std::string name;
		std::vector<double> vecistant;
	};

	private:
		std::vector<const Fem2D::Mesh *> _vecmesh;
		std::vector<tsinfo> _vecofts;
		std::string _nameoffile;
/*! This string contains the name of data file with \\ where there's a \ in the path*/
		std::string _nameofdatafile;
//! files containing the data and the timeseries
		std::ofstream _ofdata, _ofts;

/*!This function is called frequently, so if the main program crashes the files are good
 * and you need only write "end" at the end of data file: echo end>>nameoffile.data and then the files are good
 */
		void save_header () {
			std::string s = _nameoffile;

			s.append(".dx");
			_ofts.open(s.c_str(), std::ios_base::out);

			for (int i = 0; i < _vecofts.size(); ++i) {
				_ofts << "object \"" << _vecofts[i].name << "\" class series" << std::endl;

				for (int j = 0; j < _vecofts[i].vecistant.size(); ++j) {
					_ofts << "member " << j << " value file \"" << _nameofdatafile << "\",\"" << _vecofts[i].name << "_" << j << "\" position " << _vecofts[i].vecistant[j] << std::endl;
				}

				_ofts << std::endl;
			}

			_ofts << "end" << std::endl;
			_ofts.close();
		}

	public:
		DxWriter () {std::cout << "Constructor of DxWriter" << endl;}

		void openfiles (const std::string &s) {
			_nameoffile = s;
			std::string tmp = s + ".data";
			std::cout << tmp << " ";
			_ofdata.open(tmp.c_str(), std::ios_base::out);
			_nameofdatafile = "";

			for (int i = 0; i < tmp.length(); ++i) {
				if (tmp.at(i) == '\\') {
					_nameofdatafile.append(1, '\\');
				}

				_nameofdatafile.append(1, tmp.at(i));
			}
		}

		void addmesh (const Fem2D::Mesh *mesh) {
			const Fem2D::Mesh &Th(*mesh);

			_vecmesh.push_back(mesh);
			_ofdata.flags(std::ios_base::scientific);
			_ofdata.precision(15);
			_ofdata << "object \"pos_" << _vecmesh.size() - 1 << "\" class array type float rank 1 shape 2 items " << Th.nv << " data follows" << std::endl;

			for (int k = 0; k < Th.nv; ++k) {	// Scorre tutti i vertici
				_ofdata << Th(k).x << " " << Th(k).y << endl;
			}

			_ofdata << std::endl;
			_ofdata.flags(std::ios_base::fixed);
			_ofdata << "object \"conn_" << _vecmesh.size() - 1 << "\" class array type int rank 1 shape 3 items " << Th.nt << " data follows " << endl;

			for (int i = 0; i < Th.nt; i++) {
				for (int j = 0; j < 3; j++) {
					_ofdata << Th(i, j) << " ";
				}

				_ofdata << endl;
			}

			_ofdata << "attribute \"element type\" string \"triangles\" " << std::endl;
			_ofdata << "attribute \"ref\" string \"positions\" " << std::endl << std::endl;
		}

		/*!Add a new time series, defined on the mesh*/
		void addtimeseries (const string &nameofts, const Fem2D::Mesh *mesh) {
			tsinfo ts;

			ts.name = nameofts;
			std::vector<const Fem2D::Mesh *>::const_iterator first = _vecmesh.begin(), last = _vecmesh.end();

			if (std::find(first, last, mesh) == last) {
				addmesh(mesh);
				ts.imesh = _vecmesh.size() - 1;
			} else {
				ts.imesh = std::distance(first, std::find(first, last, mesh));
			}

			_vecofts.push_back(ts);
		}

		/*!Add an instant to a time series name*/
		void addistant2ts (const string &nameofts, const double t, const KN<double> &val) {
			int jj = -1;

			for (int i = 0; i < _vecofts.size(); ++i) {
				if (_vecofts[i].name == nameofts) {jj = i;}
			}

			_vecofts[jj].vecistant.push_back(t);
			_ofdata.flags(std::ios_base::scientific);
			_ofdata.precision(15);
			_ofdata << "object \"" << nameofts << "_data_" << _vecofts[jj].vecistant.size() - 1 << "\" class array type float rank 0 items " <<
			    val.size() << " data follows" << std::endl;

			for (int i = 0; i < val.size(); ++i) {
				_ofdata << val[i] << std::endl;
			}

			_ofdata << "attribute \"dep\" string \"positions\"" << std::endl << std::endl;
			_ofdata << "object \"" << nameofts << "_" << _vecofts[jj].vecistant.size() - 1 << "\" class field" << std::endl;
			_ofdata << "component \"positions\" value \"pos_" << _vecofts[jj].imesh << "\"" << std::endl;
			_ofdata << "component \"connections\" value \"conn_" << _vecofts[jj].imesh << "\"" << std::endl;
			_ofdata << "component \"data\" value \"" << nameofts << "_data_" << _vecofts[jj].vecistant.size() - 1 << "\"" << std::endl << std::endl;
			_ofdata.flush();
			save_header();
		}

		/*!Add a field*/
		void addfield (const string &nameoffield, const Fem2D::Mesh *mesh, const KN<double> &val) {
			std::vector<const Fem2D::Mesh *>::const_iterator first = _vecmesh.begin(), last = _vecmesh.end();
			int im;

			if (std::find(first, last, mesh) == last) {
				addmesh(mesh);
				im = _vecmesh.size() - 1;
			} else {
				im = std::distance(first, std::find(first, last, mesh));
			}

			_ofdata.flags(std::ios_base::scientific);
			_ofdata.precision(15);
			_ofdata << "object \"" << nameoffield << "_data\" class array type float rank 0 items " <<
			    val.size() << " data follows" << std::endl;

			for (int i = 0; i < val.size(); ++i) {
				_ofdata << val[i] << std::endl;
			}

			_ofdata << "attribute \"dep\" string \"positions\"" << std::endl << std::endl;
			_ofdata << "object \"" << nameoffield << "\" class field" << std::endl;
			_ofdata << "component \"positions\" value \"pos_" << im << "\"" << std::endl;
			_ofdata << "component \"connections\" value \"conn_" << im << "\"" << std::endl;
			_ofdata << "component \"data\" value \"" << nameoffield << "_data\"" << std::endl << std::endl;
			_ofdata.flush();
		}

		/*!Get the mesh associated with the series nameofts*/
		const Fem2D::Mesh*getmeshts (const string &nameofts) {
			for (int i = 0; i < _vecofts.size(); ++i) {
				if (_vecofts[i].name == nameofts) {return _vecmesh[_vecofts[i].imesh];}
			}

			return NULL;
		}

		void init () {
			new(this) DxWriter();
		}

		void destroy () {
			if (_ofdata.is_open()) {
				_ofdata << std::endl << "end" << std::endl;
				_ofdata.close();
			}
		}
};

class Dxwritesol_Op: public E_F0mps
{
	public:
		typedef long Result;
		Expression edx;
		Expression ename;	//!< name of time series or field
		Expression et;	//!< time
		long what;	// 1 scalar, 2 vector, 3 symtensor
		long nbfloat;	// 1 scalar, n vector (3D), n symtensor(3D)
		Expression evct;

	public:
		Dxwritesol_Op (const basicAC_F0 &args):  what(0), nbfloat(0) {
			evct = 0;
			// int nbofsol;
			// int ddim = 2;
			// int stsize = 3;
			// There's no named parameter
			args.SetNameParam();
			if (args.size() != 4) {
				CompileError("Dxwritesol accepts only 4 parameters");
			}

			if (BCastTo<DxWriter *>(args[0])) {edx = CastTo<DxWriter *>(args[0]);}

			if (BCastTo<string *>(args[1])) {ename = CastTo<string *>(args[1]);}

			if (BCastTo<double>(args[2])) {et = CastTo<double>(args[2]);}

			if (args[3].left() == atype<double>()) {
				what = 1;
				nbfloat = 1;
				evct = to<double>(args[3]);
			} else if (args[3].left() == atype<double *>()) {
				what = 1;
				nbfloat = 1;
				evct = to<double>(args[3]);
			} else if (BCastTo<pfer>(args[3])) {
				what = 1;
				nbfloat = 1;
				evct = to<double>(args[3]);
			} else if (args[3].left() == atype<E_Array>()) {
				CompileError("Until now only scalar solution");

				/*const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
				 * //cout << "taille" << a0->size() << endl;
				 * //if (a0->size() != ddim || a0->size() != stsize)
				 * //  CompileError("savesol in 2D: vector solution is 2 composant, symmetric solution is 3 composant");
				 * if( a0->size() == ddim){
				 *      // vector solution
				 *      what=2;
				 *      nbfloat=a0->size();
				 *      for(int j=0; j<l[i].nbfloat; j++){
				 *              //evct[j] = to<double>( (*a0)[j]);
				 *      }
				 * }
				 * else if( a0->size() == stsize){
				 *      // symmetric tensor solution
				 *      what=3;
				 *      nbfloat=a0->size();
				 *      for(int j=0; j<l[i].nbfloat; j++){
				 *              //evct[j] = to<double>( (*a0)[j]);
				 *      }
				 * }*/
			} else {
				CompileError("savesol in 2D: Sorry no way to save this kind of data");
			}
		}

		static ArrayOfaType typeargs () {return ArrayOfaType(atype<DxWriter *>(), atype<string *>(), atype<double>(), true);}	// all type

		static E_F0*f (const basicAC_F0 &args) {return new Dxwritesol_Op(args);}

		AnyType operator () (Stack stack)  const;
};

AnyType Dxwritesol_Op::operator () (Stack stack)  const {
	MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
	DxWriter &dx = *(GetAny<DxWriter *>((*edx)(stack)));
	string &name = *(GetAny<string *>((*ename)(stack)));
	double t = GetAny<double>((*et)(stack));
	const Mesh &Th = *(dx.getmeshts(name));
	int nt = Th.nt;
	int nv = Th.nv;
	int nbsol = nv;
	long longdefault = 0;

	KN<double> valsol(nbsol);
	valsol = 0.;
	KN<int> takemesh(nbsol);
	takemesh = 0;
	MeshPoint *mp3(MeshPointStack(stack));

	for (int it = 0; it < nt; it++) {
		for (int iv = 0; iv < 3; iv++) {
			int i = Th(it, iv);
			mp3->setP(&Th, it, iv);
			valsol[i] = valsol[i] + GetAny<double>((*evct)(stack));
			++takemesh[i];
		}
	}

	for (int i = 0; i < nbsol; i++) {
		valsol[i] /= takemesh[i];
	}

	// Writes valsol on the file file
	dx.addistant2ts(name, t, valsol);

	return longdefault;
}

// le vrai constructeur est la
DxWriter*init_DxWriter (DxWriter *const &a, string *const &s) {
	std::cout << "start init_DxWriter" << std::endl;

	a->init();
	a->openfiles(*s);
	std::cout << "end init_DxWriter" << std::endl;
	return a;
}

void*call_addmesh (DxWriter *const &mt, const Fem2D::Mesh *const &pTh) {
	mt->addmesh(pTh);
	return NULL;
}

void*call_addtimeseries (DxWriter *const &mt, string *const &name, const Fem2D::Mesh *const &pTh) {
	mt->addtimeseries(*name, pTh);
	return NULL;
}

// Add the function name to the freefem++ table
/*  class Init { public:
 * Init();
 * };
 * $1 */
static void Load_Init () {
	Dcl_Type<DxWriter *>(InitP<DxWriter>, Destroy<DxWriter> );	// declare deux nouveau type pour freefem++  un pointeur et

	zzzfff->Add("DxWriter", atype<DxWriter *>());	// ajoute le type myType a freefem++
	// constructeur  d'un type myType  dans freefem
	TheOperators->Add("<-", new OneOperator2_<DxWriter *, DxWriter *, string *>(&init_DxWriter));

	Global.Add("Dxaddmesh", "(", new OneOperator2_<void *, DxWriter *, const Fem2D::Mesh *>(call_addmesh));
	Global.Add("Dxaddtimeseries", "(", new OneOperator3_<void *, DxWriter *, std::string *, const Fem2D::Mesh *>(call_addtimeseries));

	Global.Add("Dxaddsol2ts", "(", new OneOperatorCode<Dxwritesol_Op> );

	// atype< myType * >()->Add("(","",new OneOperator3_<myType_uv,myType *,double,double  >(set_myType_uv));
}

LOADFUNC(Load_Init)
